data(faithful)
waiting = as.matrix(faithful[, 2, drop = FALSE])
data(faithful)
waiting = as.matrix(faithful[, 2, drop = FALSE])
waiting_plot <- ggplot(as.data.frame(waiting), aes(x = waiting)) +
geom_histogram(aes(y = after_stat(density)),
breaks = seq(min(waiting), max(waiting), length.out = 21),
col = "lightcyan4", fill = "lightcyan3") +
geom_density(col = "indianred2", linewidth = 1
) +
labs(x = "Time between two eruptions (in minutes)", y = "Density") +
theme_minimal() +
ggtitle("Density and Histogram of Elapsing Time Between Eruption")
waiting_plot
# prepare the data
X <- waiting
N <- length(X)
# define log.lik helper function
log.lik <- function(X, mu_1, mu_2, var_1, var_2, pi_1, pi_2) {
sum(log(pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(X, mu_2, sd = sqrt(var_2))))
}
# em function
em_mixture <- function(starting_values, X, tol = .0001, maxits = 100) {
# initialize convergence to false and iteration number to 0
converged <- FALSE
iter <- 0
N <- length(X)
# initialize starting values
pi_1 <- starting_values$pi_1
pi_2 <- starting_values$pi_2
mu_1 <- starting_values$mu_1
mu_2 <- starting_values$mu_2
var_1 <- starting_values$var_1
var_2 <- starting_values$var_2
# to save vector of log likelihoods and means
log_likelihoods <- numeric(maxits)
saved_mu_1 <- numeric(maxits)
saved_mu_2 <- numeric(maxits)
ll_changes <- numeric(maxits - 1)
saved_var_1 <- numeric(maxits)
saved_var_2 <- numeric(maxits)
saved_pi_1 <- numeric(maxits)
saved_pi_2 <- numeric(maxits)
while ((!converged) & (iter < maxits)) {
# 1. evaluate the log likelihood at the initial parameters +
ll <- log.lik(X = X,
pi_1 = pi_1,
pi_2 = pi_2,
mu_1 = mu_1,
mu_2 = mu_2,
var_1 = var_1,
var_2 = var_2)
# 2. E-Step
gamma_1 <- pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) /
(pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(X, mu_2, sd = sqrt(var_2)))
gamma_2 <- pi_2 * dnorm(X, mu_2, sd = sqrt(var_2)) /
(pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(X, mu_2, sd = sqrt(var_2)))
# 3. M-Step
#
pi_1 <- sum(gamma_1)/N
pi_2 <- sum(gamma_2)/N
mu_1 <- sum(X * gamma_1) / sum(gamma_1)
mu_2 <- sum(X * gamma_2) / sum(gamma_2)
var_1 <- sum((X - mu_1)^2 * gamma_1) / sum(gamma_1)
var_2 <- sum((X - mu_2)^2 * gamma_2) / sum(gamma_2)
# 4. evaluate the log likelihood at the new parameter values
ll.new <- log.lik(X = X,
pi_1 = pi_1,
pi_2 = pi_2,
mu_1 = mu_1,
mu_2 = mu_2,
var_1 = var_1,
var_2 = var_2)
# store new parameters
log_likelihoods[iter + 1] <- ll.new
saved_mu_1[iter + 1] <- mu_1
saved_mu_2[iter + 1] <- mu_2
saved_var_1[iter + 1] <- var_1
saved_var_2[iter + 1] <- var_2
saved_pi_1[iter + 1] <- pi_1
saved_pi_2[iter + 1] <- pi_2
if (iter > 0) {
ll_changes[iter] <- abs(ll - ll.new) }
# 5. check convergence
if(abs(ll - ll.new) < tol) { converged <- TRUE
}
# next iteration
iter <- iter + 1
# message to keep track of progress
cat(paste0("Running iteration ", iter, ". Log likelihood changed by ", round(abs(ll - ll.new), 4), "\n"))
}
# save the parameter values
params <- data.frame(pi_1 = pi_1,
pi_2 = pi_2,
mu_1 = mu_1,
mu_2 = mu_2,
var_1 = var_1,
var_2 = var_2) %>%
pivot_longer(
cols = everything(),
names_to = c(".value", "group"), names_pattern = "(.*)_(.)"
)
final_gammas_df <- data.frame( obs = 1:length(gamma_1),
waiting_1 = round(gamma_1, 5),
waiting_2 = round(gamma_2, 5) )
return(list(
params = params,
log_likelihoods = log_likelihoods[1:iter],
ll_changes = ll_changes[1:(iter - 1)],
saved_mu_1 = saved_mu_1[1:iter],
saved_mu_2 = saved_mu_2[1:iter],
final_gammas = (head(final_gammas_df)),
saved_var_1 = saved_var_1[1:iter],
saved_var_2 = saved_var_2[1:iter],
saved_pi_1 = saved_pi_1[1:iter],
saved_pi_2 = saved_pi_2[1:iter]
))
}starting_values_1 <- list(pi_1 = .5, pi_2 = .5, mu_1 = 55, mu_2 = 80, var_1 = 1, var_2 = 1)
starting_values_2 <- list(pi_1 = .5, pi_2 = .5, mu_1 = 30, mu_2 = 60, var_1 = 1, var_2 = 1)
em1 <- em_mixture(starting_values = starting_values_1, X = waiting)## Running iteration 1. Log likelihood changed by 3842.1989
## Running iteration 2. Log likelihood changed by 0.2425
## Running iteration 3. Log likelihood changed by 0.0265
## Running iteration 4. Log likelihood changed by 0.0101
## Running iteration 5. Log likelihood changed by 0.0043
## Running iteration 6. Log likelihood changed by 0.0019
## Running iteration 7. Log likelihood changed by 0.0008
## Running iteration 8. Log likelihood changed by 0.0003
## Running iteration 9. Log likelihood changed by 0.0001
## Running iteration 10. Log likelihood changed by 0.0001em1## $params
## # A tibble: 2 × 4
## group pi mu var
## <chr> <dbl> <dbl> <dbl>
## 1 1 0.361 54.6 34.5
## 2 2 0.639 80.1 34.4
##
## $log_likelihoods
## [1] -1034.288 -1034.046 -1034.019 -1034.009 -1034.005 -1034.003 -1034.002
## [8] -1034.002 -1034.002 -1034.002
##
## $ll_changes
## [1] 0.24249928394 0.02654319958 0.01006483672 0.00430799398 0.00185608116
## [6] 0.00080090007 0.00034596748 0.00014956708 0.00006469558
##
## $saved_mu_1
## [1] 54.75000 54.73771 54.69984 54.67111 54.65179 54.63909 54.63077 54.62531
## [9] 54.62173 54.61938
##
## $saved_mu_2
## [1] 80.28488 80.18128 80.14547 80.12625 80.11417 80.10628 80.10109 80.09767
## [9] 80.09541 80.09393
##
## $final_gammas
## obs waiting waiting.1
## 1 1 0.00011 0.99989
## 2 2 0.99991 0.00009
## 3 3 0.00420 0.99580
## 4 4 0.96774 0.03226
## 5 5 0.00000 1.00000
## 6 6 0.99981 0.00019
##
## $saved_var_1
## [1] 34.40750 35.54276 35.31765 35.04204 34.84594 34.71642 34.63183 34.57657
## [9] 34.54040 34.51668
##
## $saved_var_2
## [1] 31.48280 33.29762 33.79021 34.02135 34.16164 34.25301 34.31332 34.35318
## [9] 34.37949 34.39684
##
## $saved_pi_1
## [1] 0.3676471 0.3648946 0.3634577 0.3625668 0.3619894 0.3616113 0.3613630
## [8] 0.3611997 0.3610924 0.3610218
##
## $saved_pi_2
## [1] 0.6323529 0.6351054 0.6365423 0.6374332 0.6380106 0.6383887 0.6386370
## [8] 0.6388003 0.6389076 0.6389782em2 <- em_mixture(starting_values = starting_values_2, X = waiting)## Running iteration 1. Log likelihood changed by 40474.0626
## Running iteration 2. Log likelihood changed by 1.4258
## Running iteration 3. Log likelihood changed by 1.5858
## Running iteration 4. Log likelihood changed by 1.3444
## Running iteration 5. Log likelihood changed by 0.6242
## Running iteration 6. Log likelihood changed by 0.287
## Running iteration 7. Log likelihood changed by 0.1564
## Running iteration 8. Log likelihood changed by 0.105
## Running iteration 9. Log likelihood changed by 0.084
## Running iteration 10. Log likelihood changed by 0.0753
## Running iteration 11. Log likelihood changed by 0.0736
## Running iteration 12. Log likelihood changed by 0.0785
## Running iteration 13. Log likelihood changed by 0.0917
## Running iteration 14. Log likelihood changed by 0.1174
## Running iteration 15. Log likelihood changed by 0.1616
## Running iteration 16. Log likelihood changed by 0.2321
## Running iteration 17. Log likelihood changed by 0.3354
## Running iteration 18. Log likelihood changed by 0.4769
## Running iteration 19. Log likelihood changed by 0.6618
## Running iteration 20. Log likelihood changed by 0.8938
## Running iteration 21. Log likelihood changed by 1.174
## Running iteration 22. Log likelihood changed by 1.5203
## Running iteration 23. Log likelihood changed by 1.9993
## Running iteration 24. Log likelihood changed by 2.6826
## Running iteration 25. Log likelihood changed by 3.4622
## Running iteration 26. Log likelihood changed by 4.0058
## Running iteration 27. Log likelihood changed by 4.1515
## Running iteration 28. Log likelihood changed by 4.0788
## Running iteration 29. Log likelihood changed by 3.9936
## Running iteration 30. Log likelihood changed by 3.9807
## Running iteration 31. Log likelihood changed by 4.0246
## Running iteration 32. Log likelihood changed by 3.9827
## Running iteration 33. Log likelihood changed by 3.6752
## Running iteration 34. Log likelihood changed by 3.0912
## Running iteration 35. Log likelihood changed by 2.399
## Running iteration 36. Log likelihood changed by 1.7363
## Running iteration 37. Log likelihood changed by 1.1553
## Running iteration 38. Log likelihood changed by 0.6942
## Running iteration 39. Log likelihood changed by 0.3767
## Running iteration 40. Log likelihood changed by 0.188
## Running iteration 41. Log likelihood changed by 0.0886
## Running iteration 42. Log likelihood changed by 0.0403
## Running iteration 43. Log likelihood changed by 0.018
## Running iteration 44. Log likelihood changed by 0.0079
## Running iteration 45. Log likelihood changed by 0.0035
## Running iteration 46. Log likelihood changed by 0.0015
## Running iteration 47. Log likelihood changed by 0.0007
## Running iteration 48. Log likelihood changed by 0.0003
## Running iteration 49. Log likelihood changed by 0.0001
## Running iteration 50. Log likelihood changed by 0.0001em2## $params
## # A tibble: 2 × 4
## group pi mu var
## <chr> <dbl> <dbl> <dbl>
## 1 1 0.361 54.6 34.4
## 2 2 0.639 80.1 34.5
##
## $log_likelihoods
## [1] -1095.345 -1093.919 -1092.334 -1090.989 -1090.365 -1090.078 -1089.922
## [8] -1089.817 -1089.733 -1089.657 -1089.584 -1089.505 -1089.414 -1089.296
## [15] -1089.135 -1088.902 -1088.567 -1088.090 -1087.428 -1086.535 -1085.360
## [22] -1083.840 -1081.841 -1079.158 -1075.696 -1071.690 -1067.539 -1063.460
## [29] -1059.466 -1055.486 -1051.461 -1047.478 -1043.803 -1040.712 -1038.313
## [36] -1036.577 -1035.421 -1034.727 -1034.351 -1034.163 -1034.074 -1034.034
## [43] -1034.016 -1034.008 -1034.004 -1034.003 -1034.002 -1034.002 -1034.002
## [50] -1034.002
##
## $ll_changes
## [1] 1.42580723273 1.58578318422 1.34440088607 0.62422788506 0.28697213883
## [6] 0.15638497314 0.10501440105 0.08401235063 0.07528822039 0.07361061792
## [11] 0.07847914200 0.09174982788 0.11739559141 0.16164337823 0.23208819827
## [16] 0.33544312382 0.47688142788 0.66181112993 0.89381152691 1.17400170099
## [21] 1.52029439942 1.99928550981 2.68262673657 3.46222489501 4.00583768251
## [26] 4.15148012708 4.07878500674 3.99355798650 3.98070687911 4.02461514856
## [31] 3.98265459756 3.67520501228 3.09121240806 2.39895017144 1.73627559810
## [36] 1.15530292225 0.69416182073 0.37665394614 0.18803046515 0.08860539929
## [41] 0.04029663788 0.01795439269 0.00790726233 0.00345934260 0.00150753994
## [46] 0.00065542944 0.00028454811 0.00012342167 0.00005350274
##
## $saved_mu_1
## [1] 44.20000 44.99270 45.57248 45.84219 45.95954 46.03948 46.10221 46.15491
## [9] 46.20275 46.24875 46.29481 46.34275 46.39498 46.45487 46.52708 46.61744
## [17] 46.73233 46.87780 47.05950 47.28271 47.55098 47.86586 48.23165 48.65997
## [25] 49.15371 49.67466 50.15822 50.57768 50.95876 51.34081 51.74386 52.16116
## [33] 52.57092 52.95541 53.30660 53.61766 53.88038 54.09020 54.24919 54.36463
## [41] 54.44583 54.50168 54.53952 54.56491 54.58183 54.59305 54.60048 54.60538
## [49] 54.60861 54.61075
##
## $saved_mu_2
## [1] 71.14471 71.10225 71.17822 71.32847 71.47849 71.58064 71.64511 71.68783
## [9] 71.71970 71.74654 71.77130 71.79569 71.82121 71.84957 71.88314 71.92512
## [17] 71.97950 72.05078 72.14393 72.26463 72.41922 72.61440 72.85895 73.16941
## [25] 73.57151 74.08209 74.68207 75.32268 75.95841 76.56663 77.14541 77.69499
## [33] 78.20031 78.63903 79.00268 79.29692 79.52894 79.70443 79.83091 79.91840
## [41] 79.97727 80.01628 80.04195 80.05881 80.06987 80.07714 80.08191 80.08504
## [49] 80.08711 80.08846
##
## $final_gammas
## obs waiting waiting.1
## 1 1 0.00010 0.99990
## 2 2 0.99991 0.00009
## 3 3 0.00408 0.99592
## 4 4 0.96705 0.03295
## 5 5 0.00000 1.00000
## 6 6 0.99981 0.00019
##
## $saved_var_1
## [1] 0.9600000 1.0268115 0.7063628 0.5077563 0.4906262 0.5289275
## [7] 0.5806714 0.6373783 0.6972916 0.7605596 0.8287625 0.9049071
## [13] 0.9936607 1.1017980 1.2386355 1.4159806 1.6473870 1.9477409
## [19] 2.3335815 2.8214029 3.4255237 4.1654256 5.0841950 6.2334865
## [25] 7.5842371 8.9823096 10.3210570 11.6557531 13.0816774 14.6591372
## [31] 16.4057082 18.2822749 20.2397722 22.2756492 24.3720843 26.4245451
## [37] 28.2933362 29.8826164 31.1597204 32.1386510 32.8598983 33.3745051
## [43] 33.7327431 33.9776927 34.1430870 34.2538084 34.3275022 34.3763624
## [49] 34.4086749 34.4300077
##
## $saved_var_2
## [1] 179.17015 180.23681 178.98106 176.31071 173.58818 171.73370 170.57450
## [8] 169.81798 169.26369 168.80405 168.38538 167.97704 167.55413 167.08822
## [15] 166.54147 165.86308 164.99006 163.85261 162.37669 160.48010 158.07309
## [22] 155.06235 151.32685 146.62998 140.55652 132.69332 123.04488 112.20023
## [29] 101.03600 90.25977 80.14728 70.74596 62.28422 55.14939 49.48332
## [36] 45.11197 41.80520 39.38798 37.69742 36.56008 35.81374 35.32951
## [43] 35.01611 34.81290 34.68073 34.59448 34.53807 34.50110 34.47684
## [50] 34.46091
##
## $saved_pi_1
## [1] 0.009191176 0.007858950 0.010980381 0.016927210 0.022784483 0.026763999
## [7] 0.029285957 0.030970839 0.032238956 0.033315887 0.034315485 0.035305789
## [13] 0.036346204 0.037508260 0.038889551 0.040622609 0.042873693 0.045831816
## [19] 0.049707049 0.054742493 0.061209144 0.069391696 0.079663450 0.092713280
## [25] 0.109528875 0.130494476 0.154340009 0.178849101 0.202456914 0.224752634
## [31] 0.245982961 0.266232466 0.284956194 0.301436265 0.315441921 0.327106898
## [37] 0.336544468 0.343846669 0.349227855 0.353033569 0.355648036 0.357411789
## [43] 0.358589106 0.359370286 0.359886876 0.360227834 0.360452618 0.360600711
## [49] 0.360698236 0.360762444
##
## $saved_pi_2
## [1] 0.9908088 0.9921411 0.9890196 0.9830728 0.9772155 0.9732360 0.9707140
## [8] 0.9690292 0.9677610 0.9666841 0.9656845 0.9646942 0.9636538 0.9624917
## [15] 0.9611104 0.9593774 0.9571263 0.9541682 0.9502930 0.9452575 0.9387909
## [22] 0.9306083 0.9203366 0.9072867 0.8904711 0.8695055 0.8456600 0.8211509
## [29] 0.7975431 0.7752474 0.7540170 0.7337675 0.7150438 0.6985637 0.6845581
## [36] 0.6728931 0.6634555 0.6561533 0.6507721 0.6469664 0.6443520 0.6425882
## [43] 0.6414109 0.6406297 0.6401131 0.6397722 0.6395474 0.6393993 0.6393018
## [50] 0.6392376Indeed, for EM2
# for em2
em_params_2 <- em2$params
em_log_likelihoods_2 <- em2$log_likelihoods
em_saved_mu_1_2 <- em2$saved_mu_1
em_saved_mu_2_2 <- em2$saved_mu_2
em_final_gammas_2 <- em2$final_gammas
ll_changes_2 <- em2$ll_changes
# plotting log likelihoods
em_ll_df_2 <- tibble(
obs = 1:length(em_log_likelihoods_2),
em_log_likelihoods = (em_log_likelihoods_2)
)
em_ll_plot_2 <- ggplot(em_ll_df_2, aes(x = obs, y = em_log_likelihoods)) +
geom_line(col = "indianred2", linewidth = 1) +
theme_minimal()
em_ll_plot_2
# plotting change in log likelihood at each iteration
ll_change_df_2 <- data.frame(
iteration = 1:length(ll_changes_2),
ll_change = ll_changes_2
)
ggplot(ll_change_df_2, aes(x = iteration, y = ll_change)) +
geom_line() +
theme_minimal()
While for EM1:
# for em1
em_params_1 <- em1$params
em_log_likelihoods_1 <- em1$log_likelihoods
em_saved_mu_1_1 <- em1$saved_mu_1
em_saved_mu_2_1 <- em1$saved_mu_2
em_final_gammas_1 <- em1$final_gammas
ll_changes_1 <- em1$ll_changes
# plotting log likelihoods
em_ll_df_1 <- tibble(
iteration = 1:length(em_log_likelihoods_1),
em_log_likelihoods = (em_log_likelihoods_1)
)
em_ll_plot_1 <- ggplot(em_ll_df_1, aes(x = iteration, y = em_log_likelihoods)) +
geom_line(col = "indianred2", linewidth = 1) +
theme_minimal()
em_ll_plot_1
# plotting change in log likelihood at each iteration
ll_change_df_1 <- data.frame(
iteration = 1:length(ll_changes_1),
ll_change = ll_changes_1
)
ggplot(ll_change_df_1, aes(x = iteration, y = ll_change)) +
geom_line() +
theme_minimal()
Asya’s Fun
It plots the log likelihood as a function of mu_1 at a particular iteration (evaluated at the current maximizing levels of the other parameters)
as a function of mu_1
log.lik2 <- function(data, mu_1, mu_2, var_1, var_2, pi_1, pi_2) {
sum(log(pi_1 * dnorm(data, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(data, mu_2, sd = sqrt(var_2))))
}
plot.lik.tot.mu_1 <- function(iter, em) {
x_values <- seq(0, 95, length.out = 1000)
ll <- sapply(x_values, function(x) log.lik2(data = waiting,
mu_1 = x,
mu_2 = em$saved_mu_2[iter],
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
plot(x_values, ll, type = "l")
points(x_values[which.max(ll)], max(ll), col = "red", pch = 19)
}
plot.lik.part.mu_1 <- function(iter, em) {
ll <- sapply(X = seq(40, 60, length.out = 1000),
FUN = function(x) log.lik2(data = waiting,
mu_1 = x,
mu_2 = em$saved_mu_2[iter],
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
plot(seq(40, 60, length.out = 1000), ll, type = "l")
points(seq(40, 60, length.out = 1000)[which.max(ll)], max(ll), col = "red", pch = 19)
}
plot.lik.tot.mu_1(em = em2, iter = 9)
plot.lik.part.mu_1(em = em2, iter = 9)
as a function of mu_2
plot.lik.tot_mu2 <- function(iter, em) {
x_values <- seq(0, 95, length.out = 1000)
ll <- sapply(x_values, function(x) log.lik2(data = waiting,
mu_1 = em$saved_mu_1[iter],
mu_2 = x,
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
plot(x_values, ll, type = "l")
points(x_values[which.max(ll)], max(ll), col = "red", pch = 19)
}
plot.lik.part_mu2 <- function(iter, em) {
ll <- sapply(X = seq(40, 60, length.out = 1000),
FUN = function(x) log.lik2(data = waiting,
mu_1 = em$saved_mu_1[iter],
mu_2 = x,
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
plot(seq(40, 60, length.out = 1000), ll, type = "l")
points(seq(40, 60, length.out = 1000)[which.max(ll)], max(ll), col = "red", pch = 19)
}
plot.lik.tot_mu2(em = em2, iter = 14)
plot.lik.part_mu2(em = em2, iter = 14)
Given fixed parameters
plot.lik.tot <- function(iter, em) {
x_values <- seq(0, 95, length.out = 1000)
ll <- sapply(x_values, function(x) log.lik2(data = x,
mu_1 = em$saved_mu_1[iter],
mu_2 = em$saved_mu_2[iter],
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
plot(x_values, ll, type = "l")
points(x_values[which.max(ll)], max(ll), col = "red", pch = 19)
}
plot.lik.part <- function(iter, em) {
ll <- sapply(X = seq(40, 60, length.out = 1000),
FUN = function(x) log.lik2(data = x,
mu_1 = em$saved_mu_1[iter],
mu_2 = em$saved_mu_2[iter],
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
plot(seq(40, 60, length.out = 1000), ll, type = "l")
points(seq(40, 60, length.out = 1000)[which.max(ll)], max(ll), col = "red", pch = 19)
}
plot.lik.tot(em = em2, iter = 14)
plot.lik.part(em = em2, iter = 14)
total animation EM2 fixed
# complete plot
plot.lik.tot.animation <- function(em, start_iter, end_iter) {
x_values <- seq(0, 90, length.out = 1000)
data_list <- lapply(start_iter:end_iter, function(iter) {
ll <- sapply(x_values, FUN = function(x) log.lik2(data = x,
mu_1 = em$saved_mu_1[iter],
mu_2 = em$saved_mu_2[iter],
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
data.frame(iter = iter, x = x_values, ll = ll)
})
data <- do.call(rbind, data_list)
p <- ggplot(data, aes(x = x, y = ll, group = iter)) +
geom_line() +
transition_states(states = iter, transition_length = 2, state_length = 1) +
labs(title = "Iteration: {closest_state}") +
theme_minimal()
animate(p, nframes = 100, duration = 5, width = 800, height = 600)
}
# Use the function like this
plot.lik.tot.animation(em = em2, start_iter = 0, end_iter = 50)
total animation EM1 fixed
plot.lik.tot.animation(em = em1, start_iter = 0, end_iter = 10)
total animation EM2 varying mu_1
# complete plot
plot.lik.tot.mu_1_animation <- function(em, start_iter, end_iter) {
x_values <- seq(0, 95, length.out = 1000)
data_list <- lapply(start_iter:end_iter, function(iter) {
ll <- sapply(x_values, function(x) log.lik2(data = waiting,
mu_1 = x,
mu_2 = em$saved_mu_2[iter],
var_1 = em$saved_var_1[iter],
var_2 = em$saved_var_2[iter],
pi_1 = em$saved_pi_1[iter],
pi_2 = em$saved_pi_2[iter]))
data.frame(iter = iter, x = x_values, ll = ll)
})
data <- do.call(rbind, data_list)
p <- ggplot(data, aes(x = x, y = ll, group = iter)) +
geom_line() +
transition_states(states = iter, transition_length = 2, state_length = 1) +
labs(title = "Iteration: {closest_state}") +
theme_minimal()
animate(p, duration = 5, nframes = 100, width = 800, height = 600)
}
# Use the function like this
plot.lik.tot.mu_1_animation(em = em2, start_iter = 1, end_iter = 30)